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CS1 ; A Two-Dimensional  Finite  Element  Charge-Sheet  Model 
of  a Short-Channel  MOS  Transistor 

C.  L,  Wilson  and  J.  L,  Blue* 

National  Bureau  of  Standards 
Washington,  DC  20234 

Abstract 

A two-dimensional  charge-sheet  model  for  short-channel  MOS  transis- 
tors has  been  developed.  The  unique  feature  of  the  model  is  that 
the  effect  of  channel  inversion  layer  charge  is  included  as  a non- 
linear integral  boundary  condition  on  the  two-dimensional  electro- 
static field  in  the  transistor.  The  average  inversion  layer  charge 
density  and  source-drain  current  are  obtained  directly  from  the 
model  rather  than  from  the  electron  density  or  electron  quasi-Fermi 
level.  The  model  retains  all  of  the  physical  detail  of  more  com- 
plex two-dimensional  models  such  as  sensitivity  to  source-drain 
profile  shape,  channel  profile,  and  oxide  field  shape.  This  allows 
the  model  to  represent  the  changes  in  drain  current  associated  with 
short-channel  effects  while  still  allowing  simple  comparison  with 
long-channel  models.  For  long-channel  transistors,  the  results  of 
this  model  are  identical  to  Brews'  long-channel  charge-sheet  model. 
The  accuracy  of  this  model  is  verified  by  modeling  a sequence  of 
transistors  with  channel  lengths  between  4,6  and  1,1  ym.  In  short- 
channel  transistors,  effects  previously  attributed  to  high  field 
mobility  are  explained  by  simple  two-dimensional  electrostatics. 

The  simulations  produced  using  this  model  have  been  compared  to 
experimental  measurements  on  an  array  of  n-channel  MOSFETs;  the 
model  is  in  good  agreement  for  transistors  with  channel  lengths  as 
short  as  1,1  ym.  In  this  verification  process,  the  model  repre- 
sented accurately  the  onset  of  subthreshold  current,  channel- 
length-induced  threshold  voltage  offset,  and  drain-field-induced 
output  conductance  changes. 

From  studies  of  numerical  accuracy,  we  conclude  that  the  charge- 
sheet  model  can  easily  simulate  drain  current  with  an  accuracy 
which  exceeds  that  required  for  most  applications.  To  obtain  5- 
percent  accuracy  for  drain  current,  a 146  element  mesh  is  suffi- 
cient. Refinement  of  the  146  element  mesh  to  a 455  element  mesh 
gives  a current  which  is  accurate  to  0,16  percent.  Average  comput- 
er time  for  a high  accuracy  solution  is  2.5  min  on  a DEC-20,'*’ 

The  numerical  solutions  were  obtained  using  general-purpose  soft- 
ware for  solving  elliptic  partial  differential  equations.  Problems 


* NBS  Center  for  Applied  Mathematics,  Scientific  Computing  Division, 
t Certain  commercial  equipment,  instruments,  or  materials  are  identified  in 
this  paper  in  order  to  adequately  specify  the  experimental  procedure.  Such 
identification  does  not  imply  recommendation  or  endorsement  by  the  National 
Bueau  of  Standards,  nor  does  it  imply  that  the  materials  or  equipment  iden- 
tified are  necessarily  the  best  available  for  the  purpose. 
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with  exact  solutions  have  been  solved  to  test  the  correctness  and 
accuracy  of  the  codes.  Also,  the  physics  included  in  this  model 
and  the  geometry  of  the  transistor  can  be  easily  changed.  The  fi- 
nite element  method  used  allows  refinement  of  oblique  triangles 
which  is  important  in  achieving  computational  efficiency.  The 
program  is  portable  and  has  been  run  on  a DEC-20,  a VAX  11/780,  a 
Cyber  175,  and  a Uni vac  1108.* 

Key  words;  elliptic  solvers;  finite  elements;  interactive  graph- 
ics; MOS  transistor. 

1 . Introduction 

A two-dimensional  formulation  of  a charge-sheet  model  of  an  MOS  transistor 
and  an  associated  computer  program  is  presented.  This  formulation  differs 
from  the  formulation  of  De  La  Moneda  [1]  in  that  it  includes  the  effects  of 
source-drain  channel  shape  and  arbitrary  channel  doping  profile  shape.  The 
formulation  differs  from  the  formulation  proposed  by!  Brews  [2]  in  that  strong 
two-dimensional  electrostatic  effects  are  included  by  solving  Poisson's  equa- 
tion using  adaptive  finite  element  techniques  in  the  bulk  semiconductor.  In 
this  formulation,  a new  approach  that  differs  from  both  [1]  and  [2]  is  used 
to  couple  the  one-dimensional  current  flow  to  the  two-dimensional  potential. 
The  effects  of  the  charge  in  the  inverted  channel  of  the  transistor  is  in- 
cluded as  a nonlinear  boundary  condition  at  the  interface  between  the  semi- 
conductor and  the  gate  oxide.  This  formulation  allows  the  two-dimensional 
electrostatic  field  and  drain  current  to  be  calculated  without  the  inclusion 
of  an  additional  two-dimensional  current  continuity  equation.  This  results 
in  a great  improvement  in  computational  efficiency  which  allows  finite  ele- 
ment model  calculations  to  be  performed  on  a minicomputer.  This  model  is 
accurate  in  all  regions  of  inverted  transistor  operation  and  works  well  even 
in  regions  of  operation  where  the  channel  current  produces  large  changes  in 
the  electrostatic  potential  in  the  channel.  The  model  retains  the  computa- 
tional efficiency  of  single  equation  models  [3]  without  imposing  the  computer 
requirements  of  coupled  equation  models  [4,5,6].  The  complexity  and  cost  of 
the  model  lies  between  the  long-channel  one-dimensional  model  of  Brews,  which 
runs  well  on  a desk  top  calculator,  and  the  coupled  equation  models,  which 
are  best  suited  to  large  computers.  Unlike  previous  single  equation  models 
[3],  the  two-dimensional  charge-sheet  model  works  in  all  regions  of  inverted 
transistor  operation  including  weak  inversion.  Use  of  the  finite-element 
method  also  allows  the  adaptive  mesh  procedures  [7]  previously  applied  to 
MESFETs  [8]  to  be  applied  to  MOSFETs. 

In  section  two,  the  formulation  of  the  two-dimensional  charge-sheet  model  and 
the  integral  boundary  condition  at  the  semi conductor -oxide  interface  are 
discussed.  In  this  formulation,  the  total  channel  current  is  computed  as  a 
function  of  the  two-dimensional  electrostatic  potential.  The  computation  of 
total  channel  current  allows  the  effects  of  this  mobile  charge  to  be  treated 
as  an  effective  surface  charge  and  to  be  included  in  subsequent  calculations 
of  the  electrostatic  potential  as  a nonlinear  integral  boundary  condition. 

In  section  two,  the  numerical  procedures  used  to  calculate  the  electrostatic 
potential  in  two  dimensions  with  a nonlinear  boundary  condition  are  also 
discussed.  The  partial  differential  equation  is  transformed  into  a system  of 


2 


nonlinear  equations  using  adaptive,  computer-generated,  finite  elements.  The 
nonlinear  field  equations  are  solved  using  a damped  Newton's  method.  The 
resultant  system  of  linear  equations  is  solved  either  by  direct  methods  or  by 
a multigrid  iteration. 

In  section  three,  the  specific  procedures  used  to  implement  the  model  are 
discussed.  The  structure  of  the  computer  program  is  discussed  and  the  func- 
tion of  each  of  the  principal  modules  is  outlined.  The  installation  require- 
ments of  the  program  are  discussed  as  are  the  sources  of  the  software  compo- 
nents which  were  developed  by  other  authors.  The  physical  limitations  of  the 
model  are  then  discussed.  These  simplifications  are  responsible  for  the  low 
cost  of  the  model  but  may,  in  some  cases,  lead  to  modeling  inaccuracies. 
Solutions  to  these  physical  limitations  are  outlined. 

In  section  four,  the  modeling  of  an  array  of  MOS  transistors  is  discussed. 
First,  the  numerical  accuracy  of  the  model  is  discussed.  Next,  the  relation- 
ship between  long-channel  theory,  as  represented  by  the  one-dimensional 
charge-sheet  model,  and  the  two-dimensional  model  is  discussed.  Using  the 
two-dimensional  model,  it  is  possible  to  obtain  accurate  models  of  transis- 
tors with  channel  lengths  less  than  1.0  ym.  The  transistors  used  in  this 
work  have  a nominal  channel  doping  density  of  1 x 10^^  cm“^.  This  low  chan- 
nel doping  density  increases  short-channel  effects  by  increasing  the  Debye 
length  in  the  channel  of  the  transistor.  If,  for  example,  transistors  with 
channel  doping  density  of  1 x 10^^  cm”^  had  been  used,  then  short-channel 
effects  comparable  to  those  seen  in  the  1.12-ym  transistor  would  occur  in  a 
transistor  with  0.35-ym  channel  length.  In  reference  [9],  the  onset  of 
short-channel  effects  is  calculated  for  the  subthreshold  region;  a transistor 
is  classified  as  short  channel  if  excess  subthreshold  current  is  present. 

For  the  channel  doping  and  oxide  thickness  used  in  this  study,  this  short- 
channel  limit  is  at  5.6  ym;  all  of  the  transistors  studied  are  short-channel 
transistors  from  a subthreshold  current  criterion. 

2.  Theory 

A.  Two-Dimensional  Charge-Sheet  Model 

The  general  form  of  the  static  semiconductor  device  equations  is  given  by 
reference  [10]: 


V^\j;  = -p/e 

P = + P ■ 
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where  ij;  is  the  electrostatic  potential;  p is  the  net  space  charge  density;  e 
is  the  dielectric  permittivity;  is  the  net  ionized  impurity  density;  n 
and  p are  the  electron  and  hole  densities;  nj^  is  the  intrinsic  carrier 
density;  (|)p  are  the  quasi-Fermi  levels;  pp  are  the  electron 

and  hole  mobilities;  and  R is  the  net  volume  generation-recombination  rate. 
The  coordinate  system  used  in  these  calculations  is  shown  in  figure  1 . The 
source-and-drain  regions  of  the  device  are  assumed  to  be  uniform  at  the 
source  or  drain  potential  and  are  removed  from  the  calculation.  This  im- 
proves the  stability  of  the  calculation  by  removing  two  regions  in  which 
space  charge  neutrality  is  achieved  by  subtracting  a large  mobile  carrier 
density  from  a large  doping  density.  This  geometry  was  discussed  previously 
in  reference  [10].  If  recombination  is  neglected  and  the  model  is  explicitly 
restricted  to  n-channel  transistors  by  neglecting  hole  current,  then  eqs  (1) 
to  (5)  become: 


= q(n  - N^)/e 

V exp  (q(ip  - (|)^ ) /kT ) V (|)^  = 0 


(6) 

(7) 


so  that  the  electron  current  density  is  given  by: 


J = qu  n.  exp 
n ^n  1 


kT 


V(j) 
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If  eq  (6)  is  used  in  both  the  oxide  and  semiconductor  regions  of  the  transis- 
tor, then  an  interior  constraint  must  be  imposed  at  the  oxide-semiconductor 
interface : 
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where  and  are  the  permittivities  of  the  gate  oxide  and  the  semicon- 
ductor, respectively;  the  subscripts  o+,o-  represent  the  two  sides  of  the 
interface;  and  Q is  the  total  net  interface  surface  charge  density.  The 
inclusion  of  these  dielectric  effects  would  require  modification  of  eq  (1): 


V • (eVi|;)  = - p . (10) 

An  alternate  approach,  used  here,  is  to  separate  the  transistor  into  oxide 
and  semiconductor  regions  and  use  eq  (9)  as  a matching  boundary  condition. 

Direct  solution  of  eqs  (6)  and  (7)  is  possible,  and  several  methods  exist  for 
doing  this  [4,5,6].  In  surface  channel  transistors,  with  or  without  im- 
planted channels,  the  electron  current  is  concentrated  in  a thin  sheet  near 
the  semiconductor-oxide  interface.  If  eq  (6)  is  solved  in  two  dimensions  to 
obtain  i|;(x,y),  then  this  interface-associated  charge  is  given  by: 
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i’igure  1.  The  geometry  used  in  the  charge-sheet  model. 
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and  the  total  electron  current  in  the  channel  is  given  by: 
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where  3 = q/kT,  the  weighted  mobility  is  given  by  y*  = — f 


(12) 


n d 


is  the  source-drain  voltage,  L is  the  channel  length,  and  W is  the  channel 
width.  The  derivation  of  these  equations  is  given  in  Appendix  A.  The  charge 
calculated  in  eq  (12)  is  used  in  eq  (9)  to  combine  the  competing  effects  of 
gate  surface  inversion  and  inversion  charge  screening  caused  by  current.  The 
two-dimensional  charge-sheet  model  is  given  by  eqs  (6),  (9),  and  (11)  where  Q 
in  eq  (9)  is: 


2 * ^surface  "‘y' 


(13) 


and  drain  current  is  given  by  eq  (12). 

B.  Numerical  Procedures 

An  important  feature  of  this  work  is  that  the  numerical  solutions  were  ob- 
tained using  general-purpose  software  for  solving  elliptic  partial  differen- 
tial equations  (PDEs).  Although  some  efficiency  is  sacrificed,  the  benefits 
obtained  by  using  the  general-purpose  software  far  outweigh  the  loss  in  effi- 
ciency. Problems  with  exact  solutions  have  been  solved  to  test  the  correct- 
ness and  accuracy  of  the  codes  and  to  determine  optimal  strategies  and  opti- 
mal values  for  parameters  in  the  code.  The  physics  included  in  the  model  and 
the  geometry  of  the  problem  can  be  easily  changed.  The  numerical  software 
was  inspired  by  finite-element  software  of  Bank  and  Sherman  [7]  and  retains 
most  of  their  philosophy.  Currently,  a general-purpose  package  to  solve 
coupled  nonlinear  systems  of  time -dependent  partial  differential  equations  is 
being  developed.  The  current  report  uses  an  intermediate  version  of  the 
package. 


This  version  of  the  package  solves  the  PDE 

V • [a(x,y,u) Vu) ] = f 


/ 3u  9u  \ 


(14) 


in  a region  bounded  by  straight  line  segments  and  circular  arc  segments.  On 
each  segment  of  the  boundary,  u must  obey  one  of  two  types  of  boundary  condi- 
tions : 


(1) 

g(x,y,u)  = 0 

(15) 

(2) 

3u/3n  + g(x,y,u)  = 0 . 

(16) 

This  is  a fairly  general  example  of  a PDE  arising  from  conservation  of  a 
flux  and  includes  both  Poisson's  equation  and  the  hole  and  electron  continu- 
ity equations  (in  the  multiple-PDE  version). 

At  the  heart  of  the  package  is  a module  which  solves  a linear  elliptic  PDE 
using  linear  finite  elements  on  a mesh  of  triangles.  (Boundary  "triangles" 
may  have  one  curved  side.)  The  initial  triangulation  is  an  input  to  the 
package;  succeeding  triangulations  are  calculated  adaptively  by  the  package. 

For  a given  triangulation  with  M vertices,  the  calculation  proceeds  as  fol- 
lows. The  M finite-element  basis  functions  {bj^}  are  linear  on  each  trian- 
gle; bjj,  is  1 at  vertex  m and  0 at  all  other  vertices.  The  set  of  vertices 
on  type  1 boundaries  is  denoted  by  D (for  Dirichlet).  The  solution  is 
approximated  by  a sum  of  basis  functions 

M 

u(x,y)  = \ ^ a b (x,y)  . (17) 

m m 

m=1 


The  coefficients  are  determined  by  a Galerkin  method  [11];  the  error  in  the 
PDE  is  made  orthogonal  to  each  of  the  basis  functions  except  those  in  D. 


(aVu)  + f 


) 
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m~ 1 , 2 , . . , M 
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(18) 


The  remaining  conditions  on  the  b's  are  that  the  type  1 boundary  conditions 
hold  exactly  [11]. 


9 = 0,  m in  D . 

'mm  mm-' 

Equation  (18)  is  integrated  by  parts  to  yield 
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(20) 


where  the  single  integral  is  over  the  type-2  parts  of  the  boundary  only. 
Equations  (19)  and  (20)  comprise  M nonlinear  equations  in  M unknowns,  the 
coefficients 
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The  M nonlinear  equations  are  solved  by  an  iterative  process,  a damped  New- 
ton's method.  Let  v be  the  Newton  step  and  use  a fraction  of  v as  a correc- 
tion to  u;  V has  the  same  form  as  u. 


v(x,y) 


m 


b (x,y) 
m 


(21  ) 


To  calculate  v,  we  replace  u by  u + v in  eqs  (19)  and  (20)  and  linearize  in 
V,  The  linearization  of  eq  (19)  yields 


a 

m 


3g(x  ,y  ,u  ) 
m m m 

3u 


3 = 0,  m in  D , 

m 


(22) 


and  the  linearization  of  eq  (20)  yields 


3v  3f —\l 

3x  3(3u/3y)  3y /j 

(23) 


/agb  = r 
m m 

2 

Substituting  for  u and  v from  eqs  (17)  and  (21)  and  doing  the  integrals  nu- 
merically gives  a set  of  M linear  equations  for  the  correction  coefficients 
{3jj,}»  These  linear  equations  are  sparse;  there  are  typically  only  nine 
nonzero  elements  per  row,  on  average,  and  M can  be  a few  hundred  or  a few 
thousand.  A satisfactory  way  to  solve  these  equations  is  directly,  with 
Gaussian  Elimination,  using  the  Yale  Sparse  Matrix  Package  [12]. 

Unless  u is  close  to  a solution  of  the  nonlinear  finite -element  equations, 

u + V may  be  worse  than  u.  A damped  Newton's  method  is  used:  replace  u by 

u + Xv,  where  X is  chosen  so  that  max  |r_  (u  + v) I < max  |r  (u)|.  This 

m “•  m 

avoids  divergence  of  the  iterative  process. 

In  practice,  the  initial  triangulation  may  be  too  coarse  to  give  the  desired 
accuracy  or  may  need  local  refinement  to  represent  the  solution  accurately  in 
a region  of  high  curvature.  After  the  package  has  converged  to  an  approxi- 
mate solution  on  the  initial  triangulation,  the  approximate  solution  may  be 
used  to  estimate  the  error  and  produce  a new  and  finer  triangulation  adap- 
tively. There  are  various  ways  to  approximate  the  error  in  each  triangle. 

For  the  present  work,  a simple  method  is  used.  In  the  silicon,  the  PDE  is 

= -p/e  . (24) 
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Since  ij;  is  linear  in  each  triangle,  = 0.  We  estimate  the  relative  error 

in  each  triangle  by 


triangle 

The  refined  triangulation  is  obtained  by  dividing  the  triangles  with  the 
largest  estimated  errors,  such  as  all  triangles  with  estimated  error  larger 
than  1/4  of  the  largest  error. 

For  the  refined  triangulation,  the  previously  calculated  u is  a good  guess  at 
the  new  u,  so  that  few  Newton  iterations  should  be  necessary.  The  hardest 
work  is  done  on  the  initial  triangulation,  the  one  with  the  fewest  vertices. 

If  there  are  two  or  more  triangulations,  the  linear  equations  on  the  finest 
triangulation  may  be  solved  by  a multilevel  iteration  [13,7],  saving  both 
time  and  storage  space.  However,  for  the  present  problem,  sufficient  accura- 
cy could  be  obtained  with  M of  only  a few  hundred,  so  that  a direct  solution 
using  the  Yale  Sparse  Matrix  Package  was  quite  adequate. 

Using  this  procedure,  the  solution  of  eqs  (9)  to  (11)  is  started  by  first 
obtaining  a trial  solution  for  the  linearized  form  of  eq  (10),  using  the 
solution  of  eq  (11)  obtained  from  the  one-dimensional  charge-sheet  model  with 
eq  (9)  as  a boundary  condition.  The  Jacobian  of  g required  for  the  solution 
of  eq  (22)  contains  the  Jacobian  of  eq  (11).  This  Jacobian  is  approximated 
by 


= q3N(y)  . (26) 

Once  the  surface  charge  and  approximate  local  field  term  are  obtained  from 
the  one-dimensional  model  and  a local-field  condition,  the  field  term  in  eq 
(9)  associated  with  the  oxide  is  obtained  using  a fast  Poisson  solver  [14]. 
The  trial  solution  is  then  used  to  obtain  a solution  for  eq  (11)  by  quadra- 
ture and  an  additional  fast  Poisson  solution  in  the  oxide.  This  iterative 
procedure  is  repeated  for  each  Newton  step  until  the  requested  convergence  is 
obtained  or  the  specified  maximum  iteration  count  is  reached. 

3.  Computer  Implementation 

A.  Software  Modules 


The  distributed  source  of  the  CS1  program  is  divided  into  four  files: 

CS1.F0R,  SYSTEM. FOR,  PDESUB.FOR,  and  DEVLIB.FOR.  These  four  files  are  ANSI 
66  FORTRAN  IV  and  can  be  compiled  on  any  FORTRAN  compiler  which  meets  this 
standard.  The  files  CS1 .HLP,  CS1.EXP,  and  CS1 .LNK  are  also  provided. 

CS1 .HLP  provides  a set  of  user  instructions  identical  to  Appendix  B of  this 
report.  CS1 .EXP  contains  an  example  set  of  program  input  and  output  which  is 
identical  to  Appendix  C of  this  report.  The  file  CS1 .LNK  is  a sample  set  of 
linkage  control  statements  which  have  been  used  to  link  the  program  on  the 
DEC  TOPS -20  system. 
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The  file  CS1.F0R  is  the  actual  charge-sheet  model  program.  This  file  con- 
tains the  subroutines  which  control  input  and  output,  model  initial  approxi- 
mations, the  oxide  field  solution,  and  the  calls  to  the  subroutines  which 
solve  the  partial  differential  equation.  The  starting  approximation  is  iden- 
tical to  the  charge-sheet  model  originally  developed  by  Brews  in  [2].  The 
output  of  this  long-channel  model  is  printed  at  the  beginning  of  each  CS1 
calculation.  The  oxide  solution  is  calculated  using  the  fast  methods  dis- 
cussed in  reference  [14];  these  subroutines  are  in  the  file  DEVLIB.FOR.  A 
complete  copy  of  the  Swarztrauber  and  Sweet  package  is  not  included.  Only 
those  parts  of  the  package  that  are  used  by  CS1  are  included. 

The  CS1  module  of  the  program  contains  parameters  TOL  and  IP(6)  which  are 
used  to  control  the  accuracy  of  the  solution.  The  parameter  TOL  is  used  to 
control  the  iteration  accuracy  at  each  multigrid  level.  The  parameter  IP(6) 
is  used  to  control  the  number  of  multigrid  levels.  These  parameters  are  set 
to  yield  very  accurate  solutions  in  the  version  of  the  program  which  is  dis- 
tributed. This  level  of  accuracy  greatly  exceeds  that  required  for  most 
applications,  but  is  necessary,  during  installation  and  testing,  to  verify 
that  the  program  has  been  correctly  installed.  Resetting  these  parameters  to 
lower  accuracy  by  increasing  the  iteration  tolerance  and/or  decreasing  the 
number  of  multigrid  levels  will  greatly  improve  run  times. 

The  file  SYSTEM. FOR  contains  subroutines  which  are  used  in  the  other  three 
program  files  but  which  are  system  dependent.  The  subroutines  FOPEN  and 
FCLOSE  are  used  to  open  and  close  input  and  output  files.  In  the  DEC  TOPS-20 
version  supplied  on  the  tape,  the  file  FCLOSE  is  a dummy  subroutine.  The 
subroutines  ISTIME  and  IFTIME  are  used  to  start  and  stop  clocks  which  are 
used  to  measure  the  execution  times  of  the  parts  of  the  partial  differential 
solution  process.  If  these  routines  are  not  used  and  dummy  subroutines  are 
substituted,  ISTIME  must  be  less  than,  NOT  equal  to,  IFTIME.  If  these  sub- 
routines are  dummy  routines,  then  the  timing  information  in  the  output 
summary  will  be  incorrect. 

The  subroutines  PLTUTL,  CURVEA,  and  ISCALE  contain  calls  to  vector  graphics 
subroutines.  The  calls  in  the  supplied  program  can  be  used  with  the  TCS 
software  supplied  by  Tektronix  and  can  be  used  to  drive  either  Tektronix 

iff 

4010  terminals  or  Hewlett-Packard  2647  terminals.  Other  graphics  packages 
can  easily  be  substituted  if  the  necessary  initialization,  absolute  motion, 
and  absolute  vector  drawing  primitives  are  available  in  either  software  or 
terminal  firmware.  If  these  subroutines  are  made  into  dummy  subroutines, 
graphics  output  will  not  be  available  from  the  program  except  through  the 
plot  output  file.  The  data  written  into  the  plot  output  file,  file  2,  are  a 
sequence  of  non-negative  x-  and  y-coordinate  pairs.  The  first  line  of  the 
file  contains  the  minimum  and  maximum  values  of  the  x-  and  y-coordi nates. 

Each  plot  is  terminated  by  a set  (or  sets)  of  negative  coordinate  pairs.  The 
subroutines  R1MACH,  D1MACH,  and  II  MACH  are  used  to  set  machine  constants 
which  are  used  in  the  other  program  modules  to  set  machine-dependent  quanti- 
ties. The  original  versions  of  the  functions  are  discussed  in  reference  [15] 
as  are  the  procedures  used  to  calculate  these  constants  for  a type  of  com- 
puter which  is  not  presently  included.  A list  of  the  types  of  computers 
which  have  these  constants  stored  in  comment  statements,  in  the  existing 


* See  disclaimer  on  page  1 . 
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versions  of  functions,  is  given  in  table  1.  To  convert  the  program  from  any 
of  these  computer  systems  to  another,  the  data  statements  for  the  old 
computer  are  converted  to  comments,  and  the  C's  on  the  data  statements  for 
the  new  computer  are  removed o 

The  conversion  of  the  program  from  one  computer  type  to  another  should  not 
require  modification  of  any  part  of  the  program  which  is  not  in  the  SYSTEM 
module.  The  program  was  converted  from  the  DEC-20  to  the  VAX  11/780  in  a few 
minutes  of  editing  time.  This  assumes  that  some  plotting  package  is  avail- 
able and  that  timing  subroutines  are  available  on  the  target  computer 
system. 

The  module  PDESUB.FOR  contains  the  subroutines  used  to  solve  the  partial 
differential  equations.  The  theory  of  these  subroutines  has  been  discussed 
in  section  2B  above. 

The  module  DEVLIB.FOR  contains  subroutines  used  by  the  other  three  modules. 
These  subroutines  have  been  obtained  from  several  sources.  The  error  han- 
dling software,  memory  management  software,  and  machine-independent  constant 
software  were  obtained  from  [15].  The  fast  elliptic  solver  was  obtained  from 
[14].  In  the  initial  approximation  quadrature  from  [16]  and  root  finding 
routines  from  [17]  are  used.  In  the  PDESUB  module,  the  Yale  Sparse  Matrix 
Package  is  used  [12]. 

B.  Physical  Restrictions 

Three  physical  restrictions  have  been  imposed  in  the  present  implementation 
of  the  CS1  program.  First,  the  two-dimensional  electric  field  dependence  of 
the  mobility  has  not  been  included.  Second,  the  detailed  slopes  of  the 
source-and-drain  doping  profiles  have  not  been  included  in  the  model.  Third, 
the  work  function  difference  between  the  gate  material  and  the  bulk  semicon- 
ductor is  not  included. 

The  first  limitation  is  the  most  serious  problem.  In  circuit  models  and  in 
some  two-dimensional  device  simulation  programs,  the  field  dependence  of  the 
mobility  is  empirically  included  using  an  expression  of  the  form  a/(1+b*E) 
where  a and  b are  empirically  determined  parameters  and  E is  the  magnitude  of 
the  electric  field.  Recent  work  by  Cooper  and  Nelson  [18]  suggests  that  this 
model  is  not  adequate  to  characterize  the  field  dependence  of  the  mobility. 

On  the  other  hand,  the  Cooper  and  Nelson  data  do  not  include  several  effects 
which  are  known  to  affect  the  low  field  value  of  the  mobility;  these  are  the 
effect  of  bulk  doping  and  oxide  interface  quality.  Since  no  definitive  model 
for  the  two-dimensional  dependence  of  the  mobility  is  available,  none  has 
been  included.  Future  versions  of  the  program  are  expected  to  include  such  a 
mode 1 . 

The  details  of  the  two-dimensional  source-and-drain  profiles  are  included  in 
the  model  only  as  two  elliptic  abrupt  junctions.  This  allows  the  details  of 
the  majority  carrier  quasi-Fermi  level  to  be  neglected  and  greatly  simplifies 
the  program,  making  the  charge-sheet  calculation  much  more  efficient.  On  the 
other  hand  at  high  source-drain-voltages,  this  approximation  can  lead  to 
convergence  problems  and/or  to  inaccurate  modeling  of  the  potential  near  the 
drain.  In  addition,  this  model  places  the  burden  of  specifying  lateral  and 
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Table  1.  Available  Computer  Types. 


Type 

Model 

Speed  (DEC-20=1) 

Burroughs 

1700 

CDC 

5700/6700/7700 

6000/7000 

3.3  (Cyber  175) 

Cray 

1 

Data  General 

Eclipse  S/200 

Harris 

220 

Honeywell 

600/6000 

IBM 

360/370 

Xerox 

Sigma  5/7/9 

PDP-1 0 

Ka  or  Ki 

1.0  (Kl) 

PDP-1 1 
VAX 

11/780 

0.72 

Uni vac 

1100 

0.57  (1108) 

Table  2.  Transistor  Parameters. 


Name 

Value 

Channel  Lengths 

4.65,  3.20,  1.89,  1.12  ym* 

Vertical  Junction  Depth 

1.05  ym 

Lateral  Junction  Diffusion 

o 

. 

00 

Channel  Impurity  Density 

1.0  X 10^^  cm 

Source-Drain  Impurity  Density 

1.0  X 10^0  cm 

Oxide  Thickness 

51 . 5 nm 

Channel  Width 

81  ym 

* These  lengths  were  obtained  from  scanning  electron  microscope 
measurements  of  polysilicon  length  made  by  P.  Roitman. 
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vertical  junction  depths  on  the  program  user.  Other  available  two- 
dimensional  device  simulation  programs,  such  as  MINIMOS  [4],  attempt  to  in- 
clude simple  two-dimensional  process  models  in  the  device  modeling  programs. 
These  simplified  models  do  not  agree  with  more  complex  two-dimensional  models 
[19]  or  with  experimental  measurement  of  two-dimensional  junction  shape  [20]. 
If  the  user  of  the  program  wishes  to  use  these  simplified  models,  then  the 
work  of  Kennedy  and  O'Brien  can  be  used  [21].  If  more  complicated  models 
[19]  or  experimental  data  are  available,  then  these  can  be  used. 

The  least  serious  of  the  restrictions  is  that  the  work  function  difference 
between  the  gate  material  and  the  bulk  semiconductor  is  not  included.  This 
difference,  if  it  is  known,  can  be  included  in  the  gate  bias. 

4.  Device  Model  Verification 

This  section  is  divided  into  three  parts.  In  the  first  part,  the  numerical 
accuracy  of  the  model  is  evaluated.  In  the  second  part,  the  one-  and  two- 
dimensional  models  are  compared  for  a transistor  similar  to  the  transistor 
studied  by  Brews  [2] . In  the  third  part,  experimentally  measured  data  are 
compared  to  the  one-  and  two-dimensional  models.  The  transistors  used  in 
this  experimental  study  were  taken  from  an  array  of  transistors  fabricated 
using  an  SEM-compatible  transistor  array  of  the  type  discussed  in  reference 
[22].  This  pattern  has  been  incorporated  into  test  pattern  NBS-28A  [22].  In 
this  study  the  polysilicon  gate  lengths  were  measured  using  a scanning  elec- 
tron microscope,  and  the  metallurgical  channel  lengths  were  obtained  using 
estimates  of  lateral-to-vertical  arsenic  diffusion  in  reference  [19].  The 
relevant  device  simulation  parameters  are  given  in  table  2.  The  heavily 
doped  source-and-drain  regions  are  excluded  from  the  region  of  calculation 
because  the  potential  in  the  regions  is  pinned  to  the  source  or  drain  poten- 
tial plus  the  built-in  potential,  and  the  ratio  of  channel  to  source-drain 
depletion  depth  is  more  than  100  to  1. 

The  numerical  accuracy  of  two  different  adaptive  meshing  procedures  is  dis- 
cussed. The  first  procedure,  referred  to  as  the  direct  method,  uses  adaptive 
mesh  generation  to  refine  the  mesh.  After  each  adaptation,  the  resulting 
linear  equations  are  solved  by  direct  methods.  Multigrid  iteration  is  not 
used.  The  second  method,  referred  to  as  the  multigrid  method,  uses  adaptive 
mesh  generation  to  obtain  a sequence  of  more  refined  mesh  levels.  A multi- 
level iteration  over  these  levels  is  then  used  for  the  solution  of  the  linear 
equations . 

The  problem  used  as  a test  case  for  numerical  accuracy  studies  is  a 1.12-ym 
n-channel  transistor  with  0.2  V of  applied  source-drain  voltage  and  a gate 
voltage  0.7  V above  threshold.  This  transistor,  under  these  bias  conditions, 
carries  a drain  current  of  88  pA.  The  level  2 mesh  and  a contour  plot  of  the 
electrostatic  potential  for  this  calculation  are  shown  in  figures  2a  and  2b. 
The  progress  of  the  iteration  process  for  this  simulation,  using  the  direct 
method,  is  given  in  table  3.  In  this  case,  and  all  other  convergent  cases, 
the  direct  and  multigrid  methods  yield  similar  accuracies  for  the  specified 
number  of  level  1 and  level  2 iterations.  The  direct  method  is  more  robust 
and  has  converged  in  cases  where  the  multigrid  method  fails. 
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15 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 


Table  3.  Iteration  Convergence. 


Level  1 

Level  2 

Newton 

norm 

A\{; 

s 

AJ 

s 

Newton 

norm 

All; 

AJ 

s 

36.1 

3.31 

771 

24.9 

3.78 

0.36 

14.6 

1.28 

634 

24.9 

3.70 

0.34 

7.24 

1.25 

340 

3.63 

0.809 

0.87 

7.22 

1.25 

146 

0.798 

0.679 

0.37 

2.14 

1 .24 

59.9 

0.653 

0.556 

0.31 

1.78 

1.22 

23.7 

0.548 

0.488 

0.30 

1.22 

1.19 

8.8 

0.488 

0.473 

0.29 

1.15 

1.13 

3.01 

0.462 

0.447 

0.26 

1.05 

1.03 

1.01 

0.431 

0.416 

0.23 

0.94 

0.91 

0.53 

0.398 

0.385 

0.21 

0.829 

0.803 

0.43 

0.365 

0.353 

0.18 

0.722 

0.698 

0.36 

0.333 

0.323 

0.16 

0.624 

0.602 

0.28 

0.304 

0.294 

0.14 

0.537 

0.510 

0.22 

0.276 

0.268 

0.12 

0.459 

0.442 

0.17 

0.250 

0.244 

0.11 

0.392 

0.377 

0.14 

0.227 

0.221 

0.09 

0.334 

0.322 

0.11 

0.205 

0.200 

0.08 

0.285 

0.274 

0.08 

0.186 

0.181 

0.07 

0.233 

0.242 

0.07 

0.168 

0.164 

0.06 

0.207 

0.199 

0.05 

0.152 

0.148 

0.04 

0.176 

0.169 

0.03 

0.137 

0.134 

0.03 

0.150 

0.144 

0.03 

0.124 

0.122 

0.02 

0.127 

0.123 

0.02 

0.113 

0.113 

0.02 

0.108 

0.104 

0.01 

0.104 

0.103 

0.01 

0.092 

0.089 

0.01 

0.094 

0.095 

0.076 

0.078 

0.088 

0.087 
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The  convergence  of  each  iteration  is  evaluated  by  plotting  the  average  Newton 
step  size  (Newton  norm)  as  a function  of  the  iteration  number.  These  data 
are  shown  for  the  direct  case  in  figure  3.  The  corresponding  multigrid  case 
is  shown  in  figure  4.  In  each  case  the  convergence  sequence  is  divided  into 
two  regions.  For  the  first  few  iterations,  the  Newton  norm  falls  rapidly 
from  a very  inaccurate  initial  value  to  a value  of  approximately  1 , in 
normalized  units,  or  approximately  0.026  V.  The  convergence  then  proceeds 
linearly,  at  a slower  uniform  rate,  until  the  convergence  criterion  is 
reached  or  the  maximum  allowed  iteration  count  is  exceeded.  As  the  number  of 
mesh  points  is  increased,  by  increasing  the  level  of  the  mesh,  the 
convergence  rate  decreases ; it  takes  more  iterations  to  reduce  the  Newton 
norm  to  some  specified  value.  The  reduced  rate  of  convergence  is  the  result 
of  the  approximate  Jacobian  in  eg  (26)  used  in  eq  (22).  The  smallest  Newton 
norm  is  obtained  with  the  level  2 mesh  after  23  iterations.  Solutions 
obtained  with  other  channel  lengths  and  bias  voltages  using  different  numbers 
of  quadrature  points  for  the  solution  of  eq  (11)  and  different  oxide  mesh 
densities  show,  in  all  cases,  that  the  optimal  solution  is  obtained  using  the 
level  2 mesh. 

Although  the  error  evaluation  technique  discussed  above  provides  an  accurate 
measure  of  the  convergence  of  the  numerical  method,  two  other  error  evalua- 
tion techniques  are  also  useful.  In  one  method,  the  maximum  change  in  the 
surface  potential  at  the  semiconductor-oxide  interface  after  each  Newton  step 
is  compared  to  the  Newton  norm  of  that  step.  A graph  of  this  relationship  is 
shown  in  figure  5.  In  the  initial  iterations  of  each  convergence  sequence, 
the  Newton  norm  falls  rapidly,  while  only  mimimal  improvement  in  the  surface 
potential  is  made.  As  the  Newton  norm  approaches  1.0,  the  values  of  Aij^g 
and  the  norm  become  identical.  Again,  the  iteration  is  divided  into  two 
parts.  In  the  first  part,  the  two-dimensional  potential  is  refined,  and  the 
surface  potential  is  only  crudely  calculated.  This  process  converges 
rapidly.  The  second  part  of  the  process  is  dominated  by  successive  refine- 
ments of  and  the  two-dimensional  potential  is  dominated  by  surface 

potential  boundary  effects.  The  level  3 iteration  divergence  occurs  because 
the  Jacobian  in  eq  (26)  is  not  sufficiently  exact  to  allow  the  Newton  method 
to  converge  for  the  multigrid  solution  of  the  linear  equations. 

The  most  important  error  estimation  method  compares  Aipg  to  AJ^j/  the  dif- 
ference in  the  drain  current  between  the  value  at  the  present  iteration  and 
the  final  value  of  the  drain  current.  These  data  are  plotted  for  the  test 
problem  in  figures  6a  and  6b  for  both  level  1 and  level  2.  In  the  level  1 
case,  the  value  of  AJ^j  falls  slowly  as  Aij^g  decreases  from  30  to  2.  As 
Aipg  decreases  from  2 to  1 , the  current  convergence  is  very  rapid.  Finally, 
if  A\pg  < 1,  the  current  converges  quadratically  with  Aipg. 

Only  the  quadratic  region  of  convergence  is  important  for  the  level  2 case 
shown  in  figure  6b.  The  initial  error  in  the  first  region  is  represented  by 
two  points  near  A\j<g  = 3,  and  no  structure  is  visible.  The  relationship  of 
AJ^  and  Aij;g  to  the  Newton  norm  is  important  in  that  it  allows  the  conver- 
gence of  the  equation  solution  process  to  be  monitored  by  the  Newton  norm. 

For  small  Newton  norms,  less  than  one,  the  approximation  AJ^  = constant 
(Norm) is  valid.  A request  for  a 0.1  Newton  norm  (i.e.,  0.1  kT/q)  allows 
the  error  in  the  current  to  be  reduced  to  0.16  percent  and  a 0.4  Newton  norm 
produces  a 1.6-percent  error. 
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Figure  3.  Convergence  of  the  direct  method  for  the  test  problem. 
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Figure  4.  Convergence  of  the  multigrid  method  for  the  test  problem. 
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Figure  5.  Relationship  between  the  surface  potential  change  and  the  Newton 
step  size  for  the  test  problem. 
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Figure  6.  a)  Relationship  between  surface  potential  change  (kT/q)  and  drain 
current  change  (arbitrary  unit,  = 1)  for  the  level  1 mesh. 
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Figure  6.  b)  Relationship  between  surface  potential  change  (kT/q)  and  drain 
current  change  (arbitrary  unit,  = 1)  for  the  level  2 mesh. 


From  this  discussion  of  numerical  accuracy,  one  can  conclude  that  the  accura- 
cy of  this  implementation  of  the  charge-sheet  model  can  greatly  exceed  the 
accuracy  required  for  most  applications.  To  obtain  5-percent  accuracy  for 
drain  current,  a level  1 mesh  using  10  Newton  steps  and  a relative  Newton 
step  convergence  criterion  of  2 percent  are  sufficient.  Refinement  of  the 
level  1 mesh  to  level  2 gives  a which  is  more  accurate  than  required. 

Next,  the  difference  between  the  one-  and  two-dimensional  charge-sheet  model 
in  their  long-channel  limit  is  discussed.  A transistor  similar  to  the  tran- 
sistor shown  in  figure  5 of  Brews'  paper  [2]  is  considered.  The  substrate 
doping  is  10^^  cm“^,  the  oxide  thickness  is  100  nm,  and  for  the  two- 
dimensional  case  the  polysilicon  length  is  4.6  ym,  with  0.3  ym  of  lateral 
junction  diffusion  at  each  edge,  and  0.5  ym  of  vertical  junction  depth.  The 
one-dimensional  theory  used  here  is  Brews'  eq  (11);  no  simplifications  have 
been  made.  The  results  of  these  one-  and  two-dimensional  calculations  are 
shown  in  figure  7.  Although  the  transistor  is  a short-channel  transistor  in 
the  sense  discussed  in  reference  [9] , no  threshold  voltage  offset  or  change 
in  output  conductance  in  the  saturation  region  is  calculated.  The  two  calcu- 
lations are  in  good  agreement  throughout  the  region  shown.  Similar  agreement 
has  been  obtained  for  a transistor  with  8.6-ym  polysilicon,  50-nm  gate  oxide, 
and  substrate  doping  of  5 x 10^*^  cm“^.  In  the  long-channel  limit,  the  one- 
and  two-dimensional  charge-sheet  models  are  in  good  agreement.  If  the  char- 
acteristics of  the  device  being  modeled  are  not  dominated  by  short-channel 
effects  in  the  triode  and  saturation  regions,  the  one-dimensional  charge- 
sheet  model  provides  an  adequate  model  of  the  transistor  well  below  the  limit 
discussed  in  reference  [9].  This  is  true  since  the  effects  studied  in  [9] 
for  subthreshold  regions  are  not  associated  with  large  current  changes  in 
other  regions  of  device  operation. 

The  two-dimensional  charge-sheet  model  can  be  used  to  characterize  MOS  tran- 
sistors in  regions  where  both  threshold  voltage  offset  and  variable  output 
conductance  are  important.  In  the  triode  region  of  operation,  as  the  channel 
length  of  the  transistor  is  decreased,  the  two-dimensional  electrostatic 
fields  associated  with  the  source-and-drain  junctions  extend  significantly 
into  the  channel  region  and  cause  surface  inversion,  even  when  the  gate  po- 
tential is  substantially  below  threshold.  This  is  an  electrostatic  effect 
and  is  well  modeled  by  the  two-dimensional  charge-sheet  model.  This  region 
of  transistor  operation  is  shown  in  figure  8 for  transistors  with  4.65-,  3.2-, 
and  1.12-ym  channel  lengths.  For  the  4.65-ym  transistor,  both  one-  and  two- 
dimensional  theories  are  adequate  to  model  the  triode  region  and  the  high 
current  part  of  the  subthreshold  region  shown  in  figure  8.  For  the  3.2-ym 
transistor,  the  two-dimensional  model  is  in  good  agreement  with  the  measured 
data,  while  the  one -dimensional  model  is  diverging  from  these  data.  For  the 
1.12-ym  transistor,  the  two-dimensional  model  shows  a threshold  voltage  off- 
set of  0.060  V.  This  threshold  voltage  offset  results  in  a 26-percent  in- 
crease in  measured  current,  as  compared  to  the  one-dimensional  model,  when 
the  gate  voltage  is  0.25  V.  At  some  channel  length  below  1.0  ym,  these  two- 
dimensional  electrostatic  effects  will  become  so  large  that  the  channel  cur- 
rent will  no  longer  be  confined  to  the  surface  of  the  transistor  and  the 
charge-sheet  model  will  fail.  This  is  not  the  case  for  the  short-channel 
transistors  modeled  here. 
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Comparison  of  the  drain  current  calculated  by  the  one-  and  two-dimensional 
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As  the  source-drain  voltage  is  increased,  the  apparent  threshold  voltage 
offset  produced  by  two-dimensional  electrostatics  is  increased,  along  with 
the  size  of  the  depletion  region  around  the  drain.  This  effect  is  shown  in 
figure  9 for  a source-drain  voltage  of  0.5  V;  four  transistors  are  shown  with 
1.1 2-,  1.89-,  3.20-  and  4.65-pm  channel  lengths.  In  the  shortest  channel 
transistor,  at  0.25-V  gate  voltage  the  measured  drain  current  is  twice  the 
value  given  by  the  one-dimensional  theory.  Agreement  with  the  two- 
dimensional  model  is  within  3 percent.  For  the  1.89-ym  transistor,  the 
agreement  with  two-dimensional  theory  is  within  5 percent;  at  0.25-V  gate 
voltage,  the  measured  current  is  1.5  times  the  current  given  by  the  one- 
dimensional model.  For  the  long-channel  lengths,  the  measured  data  and  the 
two-dimensional  model  both  agree  with  long-channel  theory  to  better  than  1 
percent. 

In  the  saturation  region,  neither  of  the  charge-sheet  models  is  satisfactory 
for  the  1.12-ym  transistor.  The  two-dimensional  model  does  not  converge, 
because  field  reversal  effects  in  the  channel  occur  along  the  entire  channel 
length.  As  the  source-drain  voltage  is  increased,  the  depletion  region 
around  the  drain  expands  under  the  channel.  This  expanded  drain  depletion 
region  causes  current  flow  away  from  the  surface  channel  and,  at  some  voltage 
below  punch- through,  causes  the  field  in  the  channel  to  reverse.  This  field 
reversal  effect  shortens  the  effective  channel  length  and  draws  charge  away 
from  the  semiconductor-oxide  surface  charge  sheet.  The  effect  of  this  dis- 
tributed charge  is  not  adequately  represented  in  the  charge-sheet  model  and 
leads  to  convergence  problems  in  the  saturation  region  of  the  transistor. 

This  effect  decreases  with  increasing  channel  doping  and  would  not  occur 
until  0.5  pm  if  channel  doping  were  1 x 10^^  cm~^* 

Although  some  convergence  problems  occur  due  to  field  reversal  in  the 
channel,  it  is  possible  to  make  useful  comparisons  between  measured  data  and 
the  two-dimensional  charge-sheet  model  for  transistors  with  channel  lengths 
longer  than  1.12  pm.  The  transition  between  output  characteristics  which 
look  like  long-channel  characteristics  in  the  saturation  region  and  those 
which  look  like  short-channel  characteristics  occurs  between  3. 20- pm  and 

1.89- pn  channel  lengths.  The  model  results  are  compared  for  these  two  cases 
in  figures  10  and  11.  In  figure  11a  the  1.89-pm  transistor  model  is  compared 
to  measured  data.  In  figure  11b  the  experimental  data  are  compared  to  the 
long-channel  theory  of  the  one -dimensional  charge-sheet  model.  There  is  no 
similarity  between  the  long-channel  theory  and  the  measured  data.  The  long- 
channel  theory  fails  because  the  effect  of  source-drain  field-induced  channel 
inversion  is  neglected.  This  channel  inversion  is  dominated  by  the  source- 
drain  field,  since  the  calculated  threshold  voltage  offset  for  this  transis- 
tor, shown  in  figure  8,  is  not  large  enough  to  explain  this  change  in  drain 
current.  The  source-drain  component  of  the  field  produces  such  an  offset,  as 
is  shown  in  figure  9.  The  additional  source-drain-induced  offset  dominates 
the  effect  shown  in  figure  11. 

Two  interesting  physical  effects  are  also  observed  in  the  calculation  of  the 

1.89- pm  channel  transistor.  In  the  saturation  region  of  the  transistor,  a 
significant  part  of  the  total  channel  is  in  the  field  reversal  region.  If 
the  charge  associated  with  this  field  reversal  region  is  removed  from  the 
calculation  of  eqs  (11)  and  (12),  good  agreement  between  the  model  and  the 
experimental  data  is  obtained  as  shown  in  figure  11a.  This  observation  also 
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Figure  9.  Comparison  of  the  triode  region  drain  current  measured  and  calculated 
using  the  one-  and  two-dimensional  models. 


MEASURED 


28 


Figure  10.  Comparison  of  the  output  characteristics  measured  and  calculated  using 
the  two-dimensional  model  for  the  3.20-iim  transistor. 
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Figure  11.  a)  Comparison  of  the  output  characteristics  measured  and  calcu- 
lated using  the  two-dimensional  model  for  a 1.89-ym  transistor. 
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Figure  11.  b)  Comparison  of  the  output  characteristics  measured  and  cal- 
culated using  the  one-dimensional  model  for  a 1.89-ym  transistor. 


30 


explains  the  modulation  of  the  output  conductance  in  figure  11a  without  the 
inclusion  of  field  dependent  mobility.  The  effective  decrease  in  conductiv- 
ity occurs  because  the  part  of  the  channel  which  is  not  experiencing  field 
inversion  is  constrained  by  strong  electrostatic  effects  originating  in  the 
drain.  In  the  transition  region  between  the  triode  and  saturation  regions  in 
figure  11b,  the  long-channel  theory,  with  the  same  mobility,  overestimates 
the  drain  current.  Strong  two-dimensional  electrostatic  effects  at  low 
source-drain  voltages  slow  the  process  of  channel  inversion  and  at  higher 
source-drain  voltages  enhance  it,  shortening  the  channel  by  field  inversion. 

5.  Conclusions 

A two-dimensional  charge-sheet  model  of  a surface-channel  MOS  transistor  has 
been  developed.  When  this  model  is  implemented  using  adaptive  finite  element 
methods,  an  efficient  high-accuracy  transistor  simulation  program  is  possi- 
ble. The  simulations  produced  using  this  model  have  been  compared  to  experi- 
mental measurements;  the  model  is  in  good  agreement  for  transistors  with 
channel  lengths  as  short  as  1.12  ym.  In  this  verification  process  the  model 
represented  accurately  the  onset  of  sub threshold  current,  channel -length- 
induced  threshold  voltage  offset,  and  drain-field-induced  output  conductance 
changes.  These  drain-field-induced  output  conductance  changes  are  success- 
fully modeled  by  constant  channel  mobility  and  do  not  require  use  of  high 
field  mobility  effects. 
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Appendix  A 

Channel  Current  Calculation 

Following  the  general  method  presented  by  Brews  [2],  a vertical  average 
(along  the  x-direction)  of  the  electron  quasi-Fermi  level  is  calculated  over 
the  channel.  The  channel  charge  is  then  given  by: 


N(y) 


z 


(A-1 ) 


This  follows  from  eq  (2)  and  the  arguments  presented  in  reference  [2]. 
The  integral  form  of  eq  (7)  then  yields; 


V • J = 0 


which  can  be  rewritten  as: 


y* 


d(j)^(y)  I 

dy  J 


(A-2) 


(A-3) 


This  equation  has  the  solution 


z 


(A-4) 


Since  the  total  change  in  the  quasi-Fermi  level  must  equal  the  source-drain 
voltage  drop,  for  potentials  relative  to  the  source: 


z = 1 
o 


(A-5) 


Substitution  of  these  factors  into  (A-1)  yields  eq  (11).  The  current  is  then 
given  by: 


(A-6) 


which  reduces  to  eq  (12). 
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Appendix  B 
File  CS1 .HLP 


. Interactive  Execution 

The  program  is  executed  by  entering: 

RUN  CS1 .EXE 

A prompt  will  appear  for  the  next  input;  enter  a file  name 
after  each  prompt,  followed  by  a carriage  return  (denoted  <cr>) 

FILENAMEKcr> 

FILENAME2<cr> 

FILENAMES <cr> 

The  first  and  third  files  must  have  been  created  prior  to  this 
execution.  They  contain  input  data  used  by  the  program.  The 
second  file  may  be  created  during  execution.  The  FILENAMES 
must  be  legal  TOPS-20  file  names;  they  are  limited  to  eight 
characters  by  an  internal  buffer  in  CS1. 

A.  FILES 

1 ) Filel  contains  data  for  each  point  to  be 
calculated.  The  file  format  is  shown  below. 

2)  File 2 is  an  empty  file  for  output  of  plotting 
data  if  the  plots  are  to  be  saved.  This  file  name 
must  be  entered  even  if  no  plotting  is  to  be  done. 

TOPS-20  will  create  the  file  if  plotting  data  is  to 

be  stored;  other  operating  systems  may  require  the  file 
to  be  created  before  the  run  begins. 

3)  Files  contains  a table  of  channel  doping  profiles. 

B.  Examples  of  FILE  FORMATS. 

Each  point  to  be  calculated  requires  5 lines  of  input. 

Data  on  a given  line  is  free-format.  Comments  may  appear 
after  all  required  data  fields.  After  each  point 
is  calculated,  the  program  will  look  for  additional  data 
until  a source  doping  of  zero  is  found  (first  line  of  input), 
after  which  the  program  will  terminate. 

EXAMPLE  OF  Filel : Data  is  in  units  given. 

Line  1 ; 

Source  doping/cm* *3 

3.0el9 

Line  2: 

Vertical  junction  depth ( urn ), Lateral  junction  depth ( urn ), Poly-Si  length(um) 
note:  Channel  length  = Poly-Si  length-2*Lateral  junction  depths 
0.55  0.84  5.00 
Line  3: 

Drain( volts ) ,Gate(volts ) , Substrate (volts ) 
note:  Source  voltage  = 0.0. 

0.15  +1.158  0.0 
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Line  4: 

Oxide  thickness (urn) /Oxide  dielectric  constant/ Interface  charge (cm** 2) 

0.057  4.0  0.0 
Line  5: 

Channel  mobi lity( cm* *2/V-sec) /Channel  width (urn) 

400.0  20.0 
Line  6; 

Line  6 may  contain  0.0  for  stopping  program/  or  may  be 

first  line  of  data  for  next  point  to  be  calculated.  Last  line  in 

data  file  must  be  zero. 

EXAMPLE  of  File 3;  Data  for  doping  density. 

Line  1 : 

Contains  integer/  Ndop/  denoting  the  length  of  the  profile  table. 

2 

Lines  2 through  Ndop+1 : 

Each  contains  an  x (urn)  and  the  doping  density  at  x (1/cm**3) 

0.0  3.e15 

100.0  3.e15 

C.  OUTPUT  DISPLAY 

All  outputs  (except  the  drain  current)  are  in  intrinsic  Debye 
lengths  (33um)  and  thermal  volts  (25mV)/  or  units  derived  from 
these  such  as,  Thermal  volts/Debye  length. 

1 . After  program  execution/  the  following  prompts  appear  on 
the  display: 

ENTER  0/1/ 2/ 3 FOR  CURRENT  AND  NEXT  POINT/  2D  PRINTING/ 

2D  PLOTTING/  CURRENT  AND  STOPPING. 

0 — print  current.  Continue  to  next  point  in  data  file/  if  any. 

1 “ print  additional  information  about  point  just  calculated. 

2 — plot  additional  information  about  point  just  calculated. 

3 — print  current  and  terminate  execution.  Do  not  continue 

to  next  data  point  in  data  file/  if  any. 

If  response  1 or  2 is  given/  further  prompts  will  be  furnished. 

2.  OPTIONS  for  printing  (response  1 ) 

The  prompt  is 
ENTER  SUMMARY /NX /NY 

SUMMARY  — for  a summary  of  the  number  of  mesh  points  and  triangles 
created/  how  much  memory  allocated/  and  how  much  time 
used/  enter  1.  If  not  desired/  enter  0. 

NX  & NY  — how  many  points  along  X and  Y are  to  be  printed  in  the 
print  matrix.  If  NX  and  NY  are  both  negative/  a 
scaled  integer  printer  plot  of  the  potential  will  result. 
If  NX  and  NY  are  both  positive/  a rectangular  grid  is 
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placed  on  the  region;  at  each  point  X,  Y,  U (the 
potential),  UX,  and  UY  (the  derivatives  of  U in  the 
X and  Y directions)  will  be  printed.  At  grid  points  outside 
the  region,  large  numbers  will  be  printed. 

3.  OPTIONS  for  plotting  (response  2) 

The  prompt  is 

1:TRIANGLES,  2:C0NT0URS,  3:SURFACE,  4:PR0FILES 
Enter  1,  2,  3,  or  4,  to  plot  on  screen. 

Enter  101,  102,  103,  or  104,  to  save  plot  data  on  File2. 

Enter  0 to  terminate  plotting. 

Each  plot  mode  will  furnish  one  or  more  further  prompts.  In  every 
screen-plotting  case,  the  final  prompt  before  drawing  the  plot  is 
Enter  g or  G to  when  ready  to  plot 
This  gives  the  user  time  to  clear  the  screen  of  to  do  any  work 
necessary  to  put  the  terminal  in  plotting  mode. 

After  each  plot  is  drawn  on  the  screen,  the  program  will  wait  for 
a carriage  return,  and  then  re-issue  the  above  prompt.  This  gives 
the  user  time  to  clear  the  screen  or  to  do  any  work  necessary  to 
change  the  terminal  out  of  plotting  mode. 

1 ; TRIANGLE  PLOTS 

The  prompt  is 
ENTER  LEVEL 

Enter  level  number  (1,  2,  ...)  to  plot  the  mesh  of  that  level. 

Enter  0 to  terminate  triangle  plotting. 

After  the  triangle  plot  has  been  drawn,  the  program  waits 
for  a carriage  return,  and  then  reissues  the  Triangle 
Plot  prompt. 

2:  CONTOUR  PLOTS 

The  first  prompt  is 
1:FUNCTI0N,  2: ABS ( GRADIENT) 

Enter  1 to  plot  the  function  (potential),  or  2 to  plot  the 
absolute  value  of  its  gradient  (absolute  value  of  electric 
field;  units  are  thermal  volts/Debye  length). 

Enter  0 to  terminate  contour  plotting. 

The  next  prompt  is 
NUMBER  OF  CONTOURS 
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Enter  the  number  of  contour  lines  to  be  drawn.  They  will  be 
drawn  equally-spaced  between  (and  including)  the  minimum  and 
maximum  values  found. 

After  the  contour  plot  has  been  drawn,  the  program  waits 
for  a carriage  return,  and  then  reissues  the  first  Contour 
Plot  prompt. 

3:  SURFACE  PLOTS 

The  first  prompt  is 
1:FUNCTI0N,  2: ABS ( GRADIENT) 

Enter  1 to  plot  the  function  (potential),  or  2 to  plot  the 
absolute  value  of  its  gradient  (absolute  value  of  electric 
field;  units  are  thermal  volts/Debye  length). 

Enter  0 to  terminate  surface  plotting. 

The  next  prompt  is 
ENTER  LEVEL 

Enter  level  number  (1,  2,  ...)  to  make  a plot  of  the  surface 
using  the  mesh  of  the  given  level. 

Enter  0 to  return  to  the  First  Surface  Plot  prompt 

The  next  prompt  is 

ENTER  VIEWING  DIRECTION — NX,NY,NZ 

Enter  three  integers.  The  surface  plot  is  a (parallel) 
projection  of  the  surface  as  seen  from  indicated  direction. 
Thus  (0,0,1)  will  give  a top  view  (like  a Triangle  Plot), 
and  (0,-l,0)  will  give  a side  view. 

After  the  surface  plot  has  been  drawn,  the  program  waits 
for  a carriage  return,  and  then  reissues  the  first  Surface 
Plot  prompt. 

4:  PROFILE  PLOTS 

The  first  prompt  is 
ENTER  NUMBER  OF  PROFILES 

Enter  an  integer,  the  number  of  profiles  to  be  drawn. 

Enter  0 to  terminate  profile  plotting. 

The  next  prompt  is 

ENTER  INDICES  OF  PROFILE  PLANE(S),  NX  AND  NY 
Enter  two  integers.  A profile  plot  is  a slice  through  the 
solution,  perpendicular  to  the  given  direction.  Thus  (0,1) 
will  give  a slice  parallel  to  the  X axis. 

After  the  profile  plot  has  been  drawn,  the  program  waits 
for  a carriage  return,  and  then  reissues  the  first  Profile 
Plot  prompt. 


38 


2.  BATCH  EXECUTION  (for  DEC-20  only) 


A.  Submit  File 

Create  a submit  file  which  contains: 

1 . RUN  CS1 .EXE 

2.  FILENAMEI 

3.  FILENAME2 

4.  FILENAMES 

5.  (Any  other  input  used  for  the  the  equivalent  interactive 

execution) . 

B.  Output  will  go  to  file  CS1 .LOG  which  can  be 

examined  with  the  editor. 


3.  SUMMARY  OF  MODEL 


A two-dimensional  charge  sheet  model  of  a surface-channel 
MOS  transistor  has  been  developed.  When  this  model  is  implemented 
using  adaptive  finite  element  methods,  an  efficient  high-accuracy 
transistor  simulation  program  is  possible.  The  simulations  produced 
using  this  model  have  been  compared  to  experimental  measurements; 
the  model  is  in  good  agreement  for  transistors  with  channel  lengths 
as  short  as  1.1 2um.  In  this  verification  process  the  model  represented 
accurately  the  onset  of  subthreshold  current,  channel-length-induced 
threshold  voltage  offset,  and  drain-field  output  conductance 
changes. 


39 


Sample  input  data  (first  file) 


Appendix  C 
File  CS1 .EXP 


1 .E20 

0.80  0.59  3.00 
0.15  1.158  0.0 
0.0570  4.0  O.OeOO 

416.0  20.0 
1 .E20 

0.80  0.59  3.00 
0.15  1.658  0.0 
0.0570  4.0  O.OeOO 

416.0  20.0 

0.0 


Sample  input  channel  profile  (third  file) 


2 

0.0  1.E15 

200.0  1.E15 


Sample  run  using  data  shown  above  (DEC  20). 

In  this  run  all  user  input  is  labled  by  (user  input). 

@run  csl .exe 

INPUT  FILE  NAME  FOR  INPUT  DATA, 

FILE  NAME  MUST  BE  <=  8 CHARACTERS. 
cs1.dat  (user  input) 

INPUT  FILE  NAME  FOR  OUTPUT  DATA  PLOTTING 
csl .pit  (user  input) 

INPUT  FILENAME  FOR  DOPING  PROFILE  DATA 


csl.dop  (user  input) 


BETA 

= 

3.868472E+01 

PS  IB 

= 

2.871 279E-01 

LB 

= 

1 .292317E-05 

COX 

= 

6.210526E-08 

BPSO 

= 

2.101 317E+01 

BPSI 

= 

2.652973E+01 

PS(0) 

= 

6.857934E-01 

NO 

= 

9.071783E+10 

PS(L) 

= 

8.104595E-01 

NL 

= 

3.404201E+10 

ID 

= 

6.762688E-06 

10 

= 

1.181 714E+02 

I.  . . 

= 

3.564658E+01 

VG 

= 

4.479690E+01 

VSD 

= 

5.802708E+00 

LCH 

= 

5.45441  IE-02 
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4 FOR  CHANNEL  PRINTING 
1 (user  input) 
enter  summary, nx,ny 


1 -6-11 

(user  input) 

3 MULTI-GRID 

LEVELS 

VERTICES 

146 

387 

TRIANGLES 

322 

982 

STORAGE  ALLOCATION 

ALLOCATED 

USED 

MAXT 

4591 

1691 

MAXV 

1722 

617 

MAXVSM 

2296 

1153 

MAXJA 

9749 

4348 

MAXU 

1099 

1099 

MAXJU 

1099 

1062 

TOTAL 

54686 

23223 

TIMING  DISTRIBUTION 

SECONDS 

PERCENT 

GRID 

0.226 

0.118 

ASMBLY 

103.892 

54.249 

MG 

84.522 

44.135 

ADAPT 

2.870 

1.499 

TOTAL 

191 .510 

100.000 

scaled  values  from  -1.111E+01  to  2.500E+01  scale  2.500E+01 


9999 

9999 

17 

-40 

-44 

-44 

9999 

9999 

12 

-42 

-44 

-44 

9999 

9999 

-2 

-43 

-44 

-44 

89 

54 

-22 

-44 

-44 

-44 

100 

13 

-37 

-44 

-44 

-44 

94 

-1 

-41 

-44 

-44 

-44 

88 

5 

-39 

-44 

-44 

-44 

63 

38 

-29 

-44 

-44 

-44 

9999 

9999 

-12 

-44 

-44 

-44 

9999 

9999 

0 

-43 

-44 

-44 

9999 

9999 

4 

-42 

-44 

-44 

ENTER  0 FOR  CURRENT  AND  NEXTPOINT 

1 FOR  2D  PRINTING 

2 FOR  2D  PLOTTING 

3 -FOR  CURRENT  AND  STOPPING 
4 FOR  CHANNEL  PRINTING 
1 (user  input) 
enter  summary, nx, ny 
0 6 11  (user  input) 


IX  Y 


U 


UX 


UY 


1 0.00000 

2 0.01839 


0.00000  1.701 41 E+38  1.701 41 E+38 

0.00000  1.70141E+38  1.70141E+38 


1 .701 41 E+38 
1 .701 41 E+38 
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3 

4 

5 

6 

7 

8 

9 

10 

1 1 

12 

1 3 

14 

1 5 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 


0.03678 

0.00000 

9.58144E-01 

-1 .24680E+03 

-9.78029E+01 

0.05517 

0.00000 

-1.041 1 1E+01 

-1 .82800E+02 

-8.38715E+01 

0.07355 

0.00000 

-1 .10993E+01 

-1 .13871E+01 

-2.47837E+00 

0.09194 

0.00000 

-1.11075E+01 

4.58229E-01 

2.03329E-01 

0.00000 

0.01253 

1 .70141E+38 

1 .70141E+38 

1 .70141E+38 

0.01839 

0.01253 

1 .70141E+38 

1 .70141E+38 

1 .70141E+38 

0.03678 

0.01253 

-7.82947E-02 

-1 .24951E+03 

-1 .99585E+02 

0.05517 

0.01253 

-1 .07115E+01 

-8.98376E+01 

-5.09973E+00 

0.07355 

0.01253 

-1 .11026E+01 

-7.23250E+00 

2.03329E-01 

0.09194 

0.01 253 

-1 .11075E+01 

2.46640E-01 

O.OOOOOE+00 

0.00000 

0.02505 

1 .70141E+38 

1 .70141E+38 

1.701 41 E+38 

0.01839 

0.02505 

1 .70141E+38 

1 .70141E+38 

1 .70141E+38 

0.03678 

0.02505 

-3.04848E+00 

-1 .04036E+03 

-3.75744E+02 

0.05517 

0.02505 

-1 .08778E+01 

-5.80910E+01 

-6.96095E+00 

0.07355 

0.02505 

-1  .11 101E+01 

2.81 150E-01 

2.59381E-01 

0.09194 

0.02505 

-1.11075E+01 

-2.32784E-02 

O.OOOOOE+00 

0.00000 

0.03758 

1 .58477E+01 

1 .79824E+03 

4.59544E+03 

0.01839 

0.03758 

9.53071E+00 

-1 .01526E+03 

-9.01933E+02 

0.03678 

0.03758 

-7.12977E+00 

-5.72754E+02 

-2.55661E+02 

0.05517 

0.03758 

-1 .10460E+01 

-5.12203E+01 

-3.48379E+01 

0.07355 

0.03758 

-1 .11071E+01 

5.61262E-03 

-5.07924E-03 

0.09194 

0.03758 

-1 .11075E+01 

-1 .92022E-02 

O.OOOOOE+00 

0.00000 

0.0501 1 

2.20992E+01 

-1 .44968E+03 

1 .27753E+02 

0.01839 

0.05011 

1 .2761 3E+00 

-8.08550E+02 

-3.37286E+02 

0.03678 

0.0501 1 

-9.83393E+00 

-2.56518E+02 

-1 .56215E+02 

0.05517 

0.0501 1 

-1 . 10851 E+01 

-1 .14023E+01 

-8.79710E+00 

0.07355 

0.0501 1 

-1.11074E+01 

5.61262E-03 

-2.44136E-01 

0.09194 

0.05011 

-1 . 1 1075E+01 

-7.34552E-03 

O.OOOOOE+00 

0.00000 

0.06264 

2.35587E+01 

-1 .60102E+03 

1. 26531 E+02 

0.01839 

0.06264 

-3.05892E-01 

-8.09359E+02 

9.27894E+01 

0.03678 

0.06264 

-1 .02042E+01 

-1 .25966E+02 

2.96531E+01 

0.05517 

0.06264 

-1 .1 1043E+01 

-2.24213E+00 

8.89146E+00 

0.07355 

0.06264 

-1 .11076E+01 

1 .89679E-01 

-1 .47744E-02 

0.09194 

0.06264 

-1 . 1 1075E+01 

-7.34552E-03 

O.OOOOOE+00 

0.00000 

0.07516 

2.49961E+01 

-1 .49096E+03 

1 .00300E+02 

0.01839 

0.07516 

3.18066E+00 

-8.78082E+02 

4.96746E+02 

0.03678 

0.07516 

-9.27189E+00 

-3.4691 5E+02 

2.35055E+02 

0.05517 

0.07516 

-1 .10720E+01 

-1 .86051E+01 

1 .51825E+01 

0.07355 

0.07516 

-1 .11074E+01 

1 .20671E-03 

2.49914E-01 

0.09194 

0.07516 

-1 .11075E+01 

-7.34552E-03 

O.OOOOOE+00 

0.00000 

0.08769 

2.22886E+01 

1 .46264E+03 

-3.13579E+03 

0.01839 

0.08769 

1 .35906E+01 

-1 .14082E+03 

1 .08186E+03 

0.03678 

0.08769 

-5.55180E+00 

-7.1 1942E+02 

3.34245E+02 

0.05517 

0.08769 

-1 .09763E+01 

-4.73173E+01 

2.89867E+01 

0.07355 

0.08769 

-1 .1 1088E+01 

1 .20671E-03 

-1 .45270E-01 

0.09194 

0.08769 

-1 .1 1075E+01 

-1 .94337E-02 

O.OOOOOE+00 

0.00000 

0.10022 

1 .70141E+38 

1 .70141E+38 

1 .701 41 E+38 

0.01839 

0.10022 

1 .70141E+38 

1 .70141E+38 

1 .701 41 E+38 

0.03678 

0.10022 

-5.58964E-01 

-1 .21393E+03 

4.47622E+02 

0.05517 

0.10022 

-1 .06938E+01 

-8.02369E+01 

3.40412E+01 

0.07355 

0.10022 

-1 .11140E+01 

-3.81814E+00 

-9.89015E-01 

0.09194 

0.10022 

-1 . 1 1075E+01 

9.71475E-02 

O.OOOOOE+00 
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0.00000 

0.11274 

1 .70141E+38 

1 .70141E+38 

1 .701 41 E+38 

56 

0.01839 

0.11274 

1 .70141E+38 

1 .701 41 E+38 

1 .701 41 E+38 

57 

0.03678 

0.11274 

3.00789E+00 

-1 .44263E+03 

2.27028E+02 

58 

0.05517 

0.11274 

-1 .03946E+01 

-1 .50872E+02 

1 .31435E+01 

59 

0.07355 

0.11274 

-1.11073E+01 

-2.80341E+00 

2.98414E+00 

60 

0.09194 

0.11274 

-1 .1 1075E+01 

1 .12634E+00 

O.OOOOOE+00 

61 

0.00000 

0.12527 

1.701 41 E+38 

1 .70141E+38 

1.70141E+38 

62 

0.01839 

0.12527 

1 .70141E+38 

1 .701 41 E+38 

1 .701 41 E+38 

63 

0.03678 

0.12527 

4.26292E+00 

-1 .42063E+03 

1 .20125E+02 

64 

0.05517 

0.12527 

-9.91306E+00 

-2.91113E+02 

1.18065E+02 

65 

0.07355 

0.12527 

-1 .10732E+01 

-4.72232E-01 

7.43968E-01 

66 

0.09194 

0.12527 

-1 . 1 1075E+01 

-1 .97903E+00 

2.9841 4E+00 

ENTER  0 FOR  CURRENT  AND  NEXTPOINT 

1 FOR  2D  PRINTING 

2 FOR  2D  PLOTTING 

3 FOR  CURRENT  AND  STOPPING 

4 FOR  CHANNEL  PRINTING 


2 (user  input) 

1;triangles,  2:contours,  3:surface,  4;profiles 
1 (user  input) 
enter  level 

1 (user  input) 

enter  g or  G when  ready  to  plot 
g (user  input) 

(See  figure  C-1  for  graphic  output) 

Lenter  g or  g when  ready  to  plot 
g (user  input) 
enter  level 

2 (user  input) 

enter  g or  G when  ready  to  plot 
g (user  input) 

(See  figure  C-2  for  graphic  output) 
enter  g or  G when  ready  to  plot 
g (user  input) 
enter  level 

3 (user  input) 

enter  g or  G when  ready  to  plot 
g (user  input) 

(See  figure  C-3  for  graphic  output) 
enter  g or  G when  ready  to  plot 
g (user  input) 
enter  level 

0 (user  input) 

l:triangles,  2:contours,  3:surface,  4;profiles 
3 (user  input) 

1:function,  2 ;abs( gradient) 

1 (user  input) 
enter  level 

1 (user  input) 

enter  viewing  direction  — nx,ny,nz 
3 -1  1 (user  input) 
enter  g or  G when  ready  to  plot 
g (user  input) 
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(See  figure  C-4  for  graphic  output) 
enter  g or  G when  ready  to  plot 
g (user  input) 

1: function,  2 cabs (gradient ) 

1 (user  input) 
enter  level 

2 (user  input) 

enter  viewing  direction  — nx,ny,nz 

3 -1  1 (user  input) 

enter  g or  G when  ready  to  plot 
g (user  input) 

(See  figure  C-5  for  graphic  output) 
enter  g or  G when  ready  to  plot 
g (user  input) 

Icfunction,  2:abs (gradient) 

1 (user  input) 
enter  level 

3 (user  input) 

enter  viewing  direction  — nx,ny,nz 
3 -1  1 (user  input) 
enter  g or  G when  ready  to  plot 
g (user  input) 

(See  figure  C-6  for  graphic  output) 
enter  g or  G when  ready  to  plot 
g (user  input) 

Icfunction,  2 cabs (gradient) 

0 (user  input) 

Ictriangles,  2ccontours,  3csurface,  4cprofiles 

2 (user  input) 

Icfunction,  2c abs (gradient ) 

1 (user  input) 
number  of  contours 
-8  (user  input) 
umin  and  umax 

-10  25  (user  input) 
ulev 

-1.000E+01  5.000E+00  O.OOOE+00  5.000E+00 

1.500E+01  2.000E+01  2.500E+01 

enter  g or  G when  ready  to  plot 
g (user  input) 

(See  figure  C-7  for  graphic  output) 
enter  g or  G when  ready  to  plot 
g (user  input) 

Icfunction,  2c abs (gradient) 

0 (user  input) 

Ictriangles,  2ccontours,  3c surface,  4cprofiles 
0 (user  input) 

ENTER  0 FOR  CURRENT  AND  NEXTPOINT 

1 FOR  2D  PRINTING 

2 FOR  2D  PLOTTING 

3 FOR  CURRENT  AND  STOPPING 

4 FOR  CHANNEL  PRINTING 
0 (user  input) 


1 .OOOE+01 
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15-2.23085E-01 


normalized  current  3.90387E+07  current  1.61832E-05 


BETA 

= 

3.868472E+01 

PSIB 

= 

2.871279E-01 

LB 

= 

1 .292317E-05 

COX 

= 

6.210526E-08 

BPSO 

= 

2.101317E+01 

BPS  I 

= 

2.810985E+01 

PS(0) 

= 

7.266396E-01 

NO 

= 

2.659051E+1 1 

PS(L) 

= 

8.694479E-01 

NL 

= 

2.012935E+11 

ID 

— 

2.564667E-05 

10 

= 

U181714E+02 

I.  .. 

= 

1 .351853E+02 

VG 

= 

6.413927E+01 

VSD 

= 

5.802708E+00 

LCH 

= 

5.45441  IE-02 

level 

1 

1 

dpsi 

-1 .421E+01 

a 

-4.181E+06 

Umax 

3.953E+01 

2 

dpsi 

1 .567E+01 

a 

-3.541E+08 

Umax 

1 .835E+01 

3 

dpsi 

-1 .018E+00 

a 

-2.516E+08 

ximax 

8.423E+00 

4 

dpsi 

-1 .028E+00 

a 

-1 .643E+08 

Umax 

5.007E+00 

5 

dpsi 

-1 .041E+00 

a 

-1 .297E+08 

Umax 

4.235E+00 

6 

dpsi 

-1 .057E+00 

a 

-1 .233E+08 

Umax 

1 . 060E+00 

7 

dpsi 

-1 .057E+00 

a 

-1 .246E+08 

Umax 

1 . 060E+00 

8 

dpsi 

-1 .057E+00 

a 

-1 .230E+08 

Umax 

1 . 060E+00 

9 

dpsi 

-1 .014E+00 
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-9.056E+07 

Umax 

2.708E-01 

9 

dpsi 

-1 .933E-01 

a 
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0.01253 
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-9.00897E+01 
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0.07355 
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1 .93562E-01 
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0.0501 1 

-1.11075E+01 
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-1 .11075E+01 

2. 30819E-03 

O.OOOOOE+00 

0.00000 

0.06264 
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0.09194 

0.07516 

-1.11075E+01 

2.30819E-03 

O.OOOOOE+00 

0.00000 

0.08769 

2.36266E+01 

6.81750E+02 

-3.01507E+03 

0.00919 

0.08769 

2.32034E+01 

-9.481 53E+02 

1 .32428E+03 

0.01 839 

0.08769 

1 .37377E+01 

-1 .14021E+03 

1.04329E+03 

0.02758 

0.08769 

2.71884E+00 

-1 .26672E+03 

9.99422E+02 

0.03678 

0.08769 

-5.52942E+00 

-7.16435E+02 

3.34261E+02 

0.04597 

0.08769 

-9.981 66E+00 

-1 .50930E+02 

7.62547E+01 

0.05517 

0.08769 

-1 .09754E+01 

-4.75933E+01 

2.90905E+01 

0.06436 

0.08769 

-1 .10790E+01 

-3.83210E+00 

4.80330E+00 

0.07355 

0.08769 

-1 .11088E+01 

1 .06331E-01 

-1 .04312E-01 

0.08275 

0.08769 

-1 . 1 1080E+01 

8.51 575E-02 

-1 .04312E-01 
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0.10022 
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1 .701 41 E+38 

0.00919 

0.10022 

1 .70141E+38 

1 .701 41 E+38 

1 .70141E+38 

0.01839 

0.10022 

1 .70141E+38 

1 .701 41 E+38 

1 .701 41 E+38 

0.02758 

0.10022 

1 .37814E+01 

-1 .98749E+03 

9.59081E+02 

0.03678 

0.10022 

-4.89539E-01 

-1 .22193E+03 

4.45889E+02 

0.04597 

0.10022 

-8.08566E+00 

-4.40321E+02 

1 .25887E+02 

0.05517 

0.10022 
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0.06436 

0.10022 
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0.11274 

1.701 41 E+38 

1 .701 41 E+38 
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1 .701 41 E+38 

1 .701 41 E+38 

0.01839 

0.11274 

1 .70141E+38 

1 .701 41 E+38 

1.701 41 E+38 

0.02758 

0.11274 

1 .97528E+01 

-2.29747E+03 

1 .03905E+02 

0.03678 

0.11274 

3.05657E+00 

-1 .43381E+03 

2.32941E+02 

0.04597 

0.11274 

-6.61 363E+00 

-6.41455E+02 

4.29263E+01 

0.05517 

0. 1 1 274 

-1 .03770E+01 

-1 .55745E+02 

1.32586E+01 

0.06436 

0.11274 

-1 .10155E+01 

-1 .44971E+01 
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0.07355 
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113 

0.01839 

0.12527 

1 .70141E+38 

1 .70141E+38 

1 .70141E+38 

114 

0.02758 

0.12527 

2.07662E+01 

-2.12272E+03 

1 .01953E+02 

115 

0.03678 

0.12527 

4.18054E+00 

-1 .44197E+03 

9.42864E+01 

116 
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0.12527 

-5.69846E+00 
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117 
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ENTER  0 

FOR  CURRENT  AND  NEXTPOINT 

1 FOR  2D  PRINTING 

2 FOR  2D  PLOTTING 

3 FOR  CURRENT  AND  STOPPING 

4 FOR  CHANNEL  PRINTING 
0 (user  input) 

normalized  current  9.00330E+07  current  3.73224E-05  1 5-2. 2461 4E-01 

CPU  time  29:18.05  Elapsed  time  2:23:34.40 
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Figure  C-1.  , Example  of  level  1 mesh. 
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Figure  C-2.  Example  of  level  2 mesh. 
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Figure  C-3.  Example  of  level  3 mesh. 
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Figure  C-4.  Example  of  level  1 projection. 
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Figure  C-5.  Example  of  level  2 projection. 
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Figure  C-6.  Example  of  level  3 projection. 
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Figure  C-7.  Example  of  contour  plot. 
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gral boundary  condition  on  the  two-dimensional  electrostatic  field  In  the  transistor.  The  average  Inver- 
sion layer  charge  density  and  source-drain  current  are  obtained  directly  from  the  model  rather  than  from 
the  electron  density  or  electron  quasI-FermI  level.  The  model  retains  all  of  the  physical  detail  of  more 
complex  two-dimensional  models  such  as  sensitivity  to  source-drain  profile  shape,  channel  profile,  and 
oxide  field  shape.  This  allows  the  model  to  represent  the  changes  In  drain  current  associated  with  short- 
channel  effects  while  still  allowing  simple  comparison  with  long-channel  models.  For  long-channel  tran- 
sistors, the  results  of  this  model  are  Identical  to  Brews'  long-channel  charge-sheet  model.  The  accuracy 
of  this  model  Is  verified  by  modeling  a sequence  of  transistors  with  channel  lengths  between  4.6  and  1.1 
pm.  In  short-channel  transistors,  effects  previously  attributed  to  high  field  mobility  are  explained  by 
simple  two-dimensional  electrostatics. 

The  simulations  produced  using  this  model  have  been  compared  to  experimental  measurements  on  an  array  of 
n-channel  MOSFETs;  the  model  Is  In  good  agreement  for  transistors  with  channel  lengths  as  short  as  1.1  pn. 
In  this  verification  process,  the  model  represented  accurately  the  onset  of  subthreshold  current,  channel- 
length-lnduced  threshold  voltage  offset,  and  draln-f leld-Induced  output  conductance  changes. 
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From  studies  of  numerical  accuracy,  we  conclude  that  the  charge-sheet  model  can  easily  simulate  drain 
current  with  an  accuracy  which  exceeds  that  required  for  most  applications.  To  obtain  5-percent  accuracy 
for  drain  current,  a 146  element  mesh  Is  sufficient.  Refinement  of  the  146  element  mesh  to  a 455  element 
mesh  gives  a current  which  Is  accurate  to  0,16  percent.  Average  computer  time  for  a high  accuracy  solu- 
tion Is  2,5  min  on  a DEC-20, 


The  numerical  solutions  were  obtained  using  general-purpose  software  for  solving  elliptic  partial  differ- 
ential equations.  Problems  with  exact  solutions  have  been  solved  to  test  the  correctness  and  accuracy  of 
the  codes.  Also,  the  physics  Included  In  this  model  and  the  geometry  of  the  transistor  can  be  easily 
changed.  The  finite  element  method  used  allows  refinement  of  oblique  triangles  which  Is  Important  In 
achieving  computational  efficiency.  The  program  Is  portable  and  has  been  run  on  a DEC-20,  a VAX  11/780,  a 
Cyber  175,  and  a Univac  1108, 
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